
#ifndef AFTERMATH_PROBABILITY_DIST_NORMAL_HPP_INCLUDED
#define AFTERMATH_PROBABILITY_DIST_NORMAL_HPP_INCLUDED

#include "../double_constants.hpp"
#include "dist_traits.hpp"

#include <cmath>
#include <string>
#include <utility>

namespace aftermath
{
    namespace probability
    {
        // Forward declaration to define traits.
        struct dist_normal;

        template <>
        struct is_continuous<dist_normal>
        {
            typedef dist_normal distribution_type;

            static const bool value = true;
        };

        struct dist_normal
        {
            typedef dist_normal type;
            typedef double result_type;

            static const std::string name;

            static const dist_normal standard;

        private:
            double m_mu, m_sigma;
            double m_cache_variance;
            double m_cache_sigma_root_two, m_cache_pdf_scale;

        public:
            dist_normal()
                : m_mu(0.0), m_sigma(1.0),
                m_cache_variance(1.0),
                m_cache_sigma_root_two(double_constants::root_two), // sqrt(2)
                m_cache_pdf_scale(double_constants::one_div_root_two_pi) // 1 / sqrt(2 pi)
            {
            }

            explicit dist_normal(double mu, double sigma = 1.0)
                : m_mu(mu), m_sigma(sigma),
                m_cache_variance(sigma * sigma),
                m_cache_sigma_root_two(sigma * double_constants::root_two), // sigma * sqrt(2)
                m_cache_pdf_scale(double_constants::one_div_root_two_pi / sigma) // 1 / (sigma * sqrt(2 pi))
            {
            }

            double expected_value() const { return this->m_mu; }
            double variance() const { return this->m_cache_variance; }
            double standard_deviation() const { return this->m_sigma; }

            double mean() const { return this->expected_value(); }
            double stddev() const { return this->standard_deviation(); }

            double cdf(double x) const { return std::erfc((this->m_mu - x) / (this->m_cache_sigma_root_two)) / 2; }
            double pdf(double x) const { return this->m_cache_pdf_scale * exp(-(x - this->m_mu) * (x - this->m_mu) / (2 * this->m_cache_variance)); }

            bool operator ==(const type& other) const
            {
                return
                    this->m_mu == other.m_mu &&
                    this->m_sigma == other.m_sigma;
            }

            bool operator !=(const type& other) const
            {
                return !this->operator==(other);
            }
        };
        
        const std::string dist_normal::name = "norm";
        
        const dist_normal dist_normal::standard = dist_normal(0.0);
    }
}

namespace std
{
    template <>
    struct hash<aftermath::probability::dist_normal>
    {
        typedef aftermath::probability::dist_normal argument_type;
        typedef std::size_t result_type;

        result_type operator()(argument_type const& x) const
        {
            return std::hash<double>()(x.expected_value()) ^ std::hash<double>()(x.standard_deviation());
        }
    };
}

#endif // AFTERMATH_PROBABILITY_DIST_NORMAL_HPP_INCLUDED
